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ABSTRACT 


The problem of accurately replicating the parameters which 
define a given system for the purposes of implementing modern 
control strategies is important. Using an Autoregressive- 
Moving Average (ARMA) representation for the unknown system, 

a model is identified by processing input/output data to esti- 
Mate the coefficients associated with the ARMA equation. Iden- 
tification of unknown system parameters using Kalman filtering 
methods was accomplished by augmenting the state vector. [In 
this thesis the Kalman filter is formulated so that parameters 
can be identified explicitly. We call this approach the 
Adaptive Kalman Identifier (AKI). 

It is shown that the Adaptive least mean square (LMS), 
Adaptive Recursive LMS and Adaptive Lattice filters are special 
Suboptimal cases of the AKI. The convergence and modeling 
properties are compared with those of the AKI by simulation 
using various types of data. 

With minor modifications, the AKI algorithm was used to 
identify the linear and non-linear ARMA models of the phase 
locked loop (PLL). A discrete PLL using a forward Euler inte- 
gration scheme was used as a source of non-linear data. The 
AKI technique appears to enable one to discern when a potential 


non-linear system enters into its non-linear mode of operation. 
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iE ROOLETLLON 


The accomplishments in the area of microprocessor tech- 
nology in the last decade have made a noticable impact in the 
area of modern control. The desire to implement modern con- 
trol theories taking advantage of these advancements, has 
made it imperative that the engineer attempt to mathemat- 
ically replicate the system he ultimately desires to control. 
Efforts in this regard have, hence, generated a growing 
interest in the area of system identification [Ref. 1,2,3,4]. 
By implication this thesis concerns itself with discrete/ 


digital signal processing. 


A. BACKGROUND 

System identification or modeling can be accomplished by 
innovative application of existing techniques which were 
generally considered filtering or state estimation methods 
(Ref. 5]. Previous research efforts have, for obvious reasons, 
focused on linear modeling; however, there iS a riSing interest 
in non-linear modeling methods [Ref. 6,7,8]. 

An attractive form for the representation of an unknown 


System is the Autoregressive-Moving Average (ARMA) equation, 


y(k) = )} a,u(k-i) - } by (k~1) (ie) 
i=0 i=l 
which states that the present output, y(k), is a linear com- 
bination of past outputs, y(k-i), and of past and present 
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inputs, u(k). Its attractiveness lies in its linear charac- 
ter and easily implementable structure uSing microprocessor/ 
computer algorithms. Its relationship to the state space 
representation of a linear system has already been established 
[Ref. 9,10] making it germane to consider that a system is 
identifiable by an equation of the form (1.1). The system 
identification problem thus entails identifying the coeffi- 
cients as and b.. 

There are several methods for computing the coefficients, 
as and b., however, it is not the intent of this brief intro- 
duction to attempt to develop even the majority of them. 
Nevertheless, it is practical to present a referenced history 
of those methods encountered in this thesis. 

Adaptive algorithms for the purpose of estimating the 
coefficients of (1.1) have always been of interest. Widrow, 
using a least means square error criterion and implementing 
the method of steepest descent, developed an adaptive algorithm 
which estimated the coefficients of the moving average process 
associated with equation (1.1) [Ref. ll]. That is, the moving 


average model, 


y(k) = agul(k) +aj,u(k-1) + a5u(k-2) ene ey 
where, 

ay = a 

ay = (a, - ab,) 


Le 





eel 


ri ayb.,) - b, - ab.) 


a, = (a eel 


il 


was identified. The LMS theory has since been extended to 
include block LMS filtering methods [Ref. 12,13,14] using 
various search techniques [Ref. 15,4(Chapter 5)]. These 
methods have enjoyed much popularity in the area of Linear 
Prediction and digital speech processing [Ref. 16,17]. 

The shortcoming of representing equation (1.1) by its 
equivalent moving average model (1.2) lies in the practical 
aspect of its implementation. That is, the infinite series 
represented by equation (1.2) must be truncated at some point 
resulting in an approximation which may not adequately repre- 
sent equation (1.1). Hence, efforts to adaptively estimate 
the coefficients for both the autoregressive (b.) and moving 
average (a. ) processes continued [Ref. 18,19,20] with various 
degrees of success. Feintuch's Adaptive Regressive LMS pro- 
cedure [Ref. 18] is still a controversial issue [Ref. 19,20]. 

From yet a different direction, Anderson and Moore suggested 
that the Kalman filtering algorithms can be used as a means 
to identify the a. and b. coefficients* Gane Cilia teLomy ell) 


Mec. 21, pp. 50-52]. This thesis exploits this application. 


anderson and Moore in fact formulate the technigue to 
compute the coefficients a,b. fe eZ, see p where 
p = nt+m. 


is, 





Pee = LNVESTIGATIONS AND CONTRIBUTIONS 

This thesis formulates an adaptive application of the Kal- 
man filter to identify the coefficients of the general ARMA 
mem@action (1.1). It is shown that the LMS adaptive filter 
[Ref. 11] and the Adaptive Recursive LMS filter [Ref. 18] 
are (1) special cases of the Adaptive ARMA Kalman identifier 
and (2) sub-optimal with respect to the underlying least 
means square (LMS) error criterion upon which the LMS adap- 
tive filter, the adaptive recursive filter and the Kalman 
filter are based. It is also demonstrated that the adaptive 
Kalman identifier has excellent convergence properties. 

By comparison, it is noted that the Kalman algorithm 
accounts for measurement noise where the LMS algorithms do 
not, a heretofore unapproached problem by LMS algorithms. 

The results indicate that the suggested modification by 
Smeereiths (Ref. 22] of the convergence factor, Keo, for the 
Meemadaptive filter is justified. 

imac application of the Kalman filter algorithm is ex- 
tended to identification of the coefficients associated with 
a special case of the generalized non-linear ARMA model [Ref. 
8]. The results of the non-linear simulations suggest a tech- 
nique for determining when a potential non-linear system 
enters its non-linear operating regime. Such a technique can 
prove valuable when on-line performance evaluation of a known 
non-linear system is required. 

Lastly, the connection between the Kalman filter algorithm 


and the lattice filter algorithm [Ref. 6] is made. An example 


16 





which demonstrates and compares the performance of both 


algorithms is given. 


C. ORGANIZATION 

Chapter II presents and exploits the Kalman filter equa- 
tions emphasizing its connections with the Yule-Walker and 
the discrete Weiner-Hopf equations. The theory is further 
developed to investigate the general ARMA case. Chapter III 
pursues the theoretical comparisons between the discrete 
Weiner-Hopf equation upon which the adaptive LMS filter is 
based and the MA form of the Kalman filter. The theoretical 
comparison between the Adaptive Recursive LMS filter and the 
general ARMA form of the Kalman filter is made in Chapter IV. 
The software methods by which linear and non-linear synthetic 
data is generated are discussed in Chapter V. A short user's 
description of the data processing programs and the options 


provided is given in Chapter VI. A non-linear application 


Geeetcne identification of the parameters of a non-linear plant 


is presented in Chapter VII followed by a discussion and pre- 


sentation of the results in Chapter VIII. 
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ce OnMOnAPRION OF THE PROBLEM 


Identifying the coefficients of the Autoregressive- 


Moving Average (ARMA) equation, 


mam) + Byy (k-1) “of re bY (k~n) = ayu(k) =F a,u(k-1) 
Sy ae a ju (k-m) (2.la) 
m n 
me USCS 2g EUR SI) ) b.y (k~i) (2a) 
i=0 i=l 


can be formulated as an adaptive Kalman identification problem. 
That is, instead of using the well known Kalman Filter equa- 
tions [Ref. 23,24] to recursively estimate the states of a 
system, one can utilize the Kalman filter equations to adap- 
tively estimate the coefficients of either an Autoregressive 
(AR), Moving Average (MA), or Autoregressive-Moving Average 
(ARMA) process by proper definition of the quantities involved. 
We call this the Adaptive Kalman Identifier for obvious 


PeeeOns . 


fae DETERMINISTIC CASE 

Consider for the moment that the a. and b. ame cOnStane, 
Then by collecting a sufficient number of measurements of the 
input u(k) and the output y(k), one can readily solve a set 
of linear equations for the a. and b. coefficients. The 


“sufficient number" that is needed is n+m+l which define the 


kis 





n+m+l linear equations which solve the n+m+l unknown coeffi- 





cients. The matrix equation to be solved takes the form: 
y (0) u(0) 0 seis 0 1 0 0 
j 
y (2) u(1) (0) ! y (0) 0 
| 
| 
j 
| 
{ 
| 
y (m) = {u(m u(m-12) u(O)}; y(m-1) y(m-2) y(0) 0 
we nnn nr nnn nnn nr nnn nr en ef ae a oe eee om eee oe eee oe eee om eee wees oe oe ee 
y (+L) u(m+l) u(l) ;y(m) —s-y@l)~—sy(1)_— 0 
| 
j 
j 
| 
| 
y (mtn) u(mtn) ... u(m) j y(n-1) y(n-2) y (0) 
a 
y = TI|--- (2.3) 
b 


It 1S interesting to note that the (n+m+l) x (n+m+l) matrix, 
ies Dlock symmetric. The solution to (2.3) is readily 


apparent, namely, 


—— = a 4 (24) 


-1 . . 
where T 1s the inverse matrix of T. Since we have assumed 
that the elements of T are perfect noiseless measurements, 


and that we somehow know a-priori the number of unknown 


Lg 





G@eerricients, then T is of full rank, namely rank (T) = 


nt+m+l, and the inverse of T exists. 


Eee lhe NONDETERMINISTIC CASE 

Rarely can the coefficients be modeled so ideally. A 
more prudent and realistic model admits that the ARMA equa- 
tion coefficients are subject to random perturbations. Further, 
it can be said that in general measurement devices introduce 
noise into the measurement data. Hence, the measurement 
model should reflect this fact. Developing the Adaptive Kalman 


Identifier along the same lines as [Ref. 21] we let, 


a, (k+1) a. (k) - w, (k) Oo ne oe TM (25a) 


een? b.{k) + Has) Sine eres) Th Pan Set) 


where the iw, (x) ,w, (k)} are samples from a zero mean, white, 
gaussian random process. Additionally, we assume that the 


noise sources are uncorrelated, 


mmye (k)w (k)}) = O for rf#7s (2.6a) 


(kK) w) Cogs O 8 tor xref s&s ( 2..'otey) 
og S 


E lw. (kK) w) (kK ye] OL tomtalileur,is (26) 


1g Ss 


Similarly, allowing for noisy meaSurement devices, the 


measurement equation (2.1b), is modelled as, 


y(k) = 
ah 


ic 


ead 


a,u(k-1) - 


Bovi(k=a)  v(k) a7) 
0 i 2 


il 





where the {v(k)} are samples from a zero mean, white, gaussian 


random process. It is also assumed that, 


Elw, (k) v(k) ] Oe cos fall 1. (2.8a) 


E lw. (k) walk] Omatorealies ( 2@8)5) 


That is, it is assumed that the noise perturbations associated 
with one coefficient are independent of the noise perturba- 
tions associated with any other coefficient, and that the 
measurement noise and the coefficient perturbations are 
independent. 

At this point it is ncessary to use judgment and experi- 
ence and utilize all the information known about the physical 


system which equation (2.7) represents to assign variances for 


the random processes tw. (k)}, Tas WS) anceeiviky }. Let 
w, (k) 
a= £E lw, (k) w. (k) ] (2.9a) 
w. (K) 
Q1 0 
0 = (2.9b) 
0 Q. 
where, Q, = di (07 ) =O (o* ) d 
; is 1ag We ; Q. = dlag Ww, and, 
me E{v(k) vik) ] (2.9c) 


and E{x] is the expected or average value of x. 


Al, 





Instead of using x's denoting the commonly used notation 
for state variables, we retain the flavor of the problem by 
using the a,'s and b.'s as the "states" of the Adaptive Kalman 
Identifier. It is hoped that a more meaningful understanding 
of the Adaptive Kalman Identifier may thus be gained. There- 


fore, define the state vector, 


ay (kK) 


a, (k) 


ey) 
= 
| 


—— = --- 2 80) 


on 
* 
[oOo 


Bk) 


Combining equations (2.5) and (2.10), we have the discrete, 


first order Gauss-Markov process [Ref. 3], 


a(k+1) a (k) 
a = S2s4 + w(k) eee 1) 
b(k+1) b (Kk) 
where, 
w. (k) 
1 
as) - ~@ | ----— : (C2) 
w. (k) 
=) 
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Pemescthatl wik) 1s an ntm+tl vector whose elements are a con- 
catenation of the noise sequences associated with the a. and 
b.. To complete the Adaptive Kalman Identifier formulation 


we define the measurement vector, 
Meeeeee—= (u(k) u({k=1) ... u({k-m) -y(k-1) ... -y(k-n) ) (27S) 


and its associated measurement equation, 


a (k) 
y(k) = H(k) Si see OS) (2a 14) 
b (k) 


Then the solution to the Adaptive Kalman Identifier problem 


Weert. 21} is, 


a(k+1|k) a(k|k-1) 
-------- = [I -K(k)H(k)] | --------| + K(k)y (k) (2.15a) 
b(k+1|k) b (k |k=1) 
T i -1 
K(k) = P(k|k(1) H (k) (H(k) P(k|k-1) H’(k) +R] (2.15b) 
P(k+1l|k) = P(k|k-1) - K(k)H(k)P(k|k-1) + Q. (2.15c) 


Equations (2.15) are initialized by assuming an initial value 
Meretene coefficients (2.10) and assigning to our assumption 
a measure of our confidence in the initial guess. That 1s, 


we pick the values, 


------- and P(0|-1) (2.16) 


Z3 





where P(0|-1) is defined as the error covariance of the 


coefficients, 


a (k) a (k |k) a(k) 
P(k|k-l) = E 352 ee ise 
b (Kk) b (k |k) b (k) 
a (k |k) i 
= Pesca (2.17) 
b (k |k) 


for k = 0. Equation (2.15b) is generally referred to as the 


Kalman gain. Two special cases are of interest: case (1); 
a = Q for all i, and case (2); b. —wOeemenr sa 1. i). 
1. Autoregressive Form (a. = 0) 


For case 1, equation (2.7) takes the form, 


mae) = - 


b,y(k-i) + v(k). (2.18) 
ti 


Il 15 


{ 


This 1s a recursive equation stating that the present output 
1S a linear combination of past outputs corrupted with addi- 
tive gaussian noise, v(k). Equation (2.18) is more formally 
recognized as an autoregressive prcoess of order n [{Ref. 25, 
26]. Writing the Kalman solution for the coefficients of the 


AR process, 


b(k+1|k) = [I -K(k)K(k)]b(k|k-1) + K(k) y(k) (2.19a) 


K(k) = P(k|k-1) H2(k) [H(k)P(k|k-1) H’ (k) +R]72 (2.19b) 
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P(k+1|/k) = P(k|k-1) - K(k)H(k)P(k|k-1) + Q. Poe) 


Alternate forms of the Kalman equations given by Maybeck 


[Ref. 27] are, 


b (k+1|k) = P(k|k-1) P(k|k-1) b (k |k-1) 
T = 
+ P(ktl/k) HT (k) Ro+y(k) (2.20a) 
P(k+1|k) = (P(k{k-1))71+H7(k) Ro? H(k)7?. (2.20b) 


We can model the fact that we have no a-priori knowledge 


about the initial values of the coefficients by letting, 


Peoleniees == 0), (om 


Further, 1£ we remain ignorant and totally doubt our previous 


estimate, then P(k|k-1) is modeled as, 


Race =e 0. (Domes 


Bouations (2.20) reduce to, 


ik 


b(k+1|k) = (H(k) R- H?(k)]7? 


= (Ge) 7 Oe (2.23) 


the weighted least squares estimate, previously encountered 
Meeting [Rer. 28,29,30], of the coefficients, b. Carry- 
ing the analysis further and letting eS = 1, we have the 


Penrose pseudo-inverse solution [Ref. 31,32], 


b(k+1|k) = [H(k) H2(k)]~+ H(k)y(k) (2.24) 
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Swerling [Ref. 28] has shown that if the weighting matrix 
in (2.23) is not the inverse of the covariance matrix of 
the measurement errors, then the accuracy of the estimated 
coefficients (2.24) will be degraded. 

The aforementioned notwithstanding, we press further 
into the analysis of equation (2.24). The product H(k)H* (x) 


can be written as, 


aca) 
SOURS”) 
H(k)H'(k) = [y(k-1) y(k-2) ... y(k-n)] (2.25) 
y (ken) 
y* (k=1) y (k-l)y (k=2) ... y(kel)y (k=n) 
y (k=-2) y (k-1) va (k=2) 
H{k)H- (k) = We 26) 


y (k-n) y (k-1) — we. y(Ken) y (k-n) 


As we let k +> ~ the expected value of (2.26) becomes the 


autocorrelation matrix, 
@iddm H(k) H>(k)} = R_(k) (2.27) 


ee yy 


where, 
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R (0) R (1) we. R 
yy °°? yy ‘*? yy ‘P) 
=) Ale 
yy (TE) Rye (0) 3 
jo | 2.28 
a ) ( ) 
Ry ony ++ Par By 60) 


Similarly, the product H(k)y(k) as k + ~ becomes 


yy oD 
Ryy (2) 

kK = DIDS 

oy ) ( ) 
Boag (Sie 


The steady state solution for the estimate of the coeffi- 
cients b(k+1l|k) is, 


-1 
k+1|k) = R iS) ae Bee0 
(k+1[k) yy 9) ona ) ( ) 


b 

—ss 
MemtietrtoOn (2.30) is one of the starting points from which 
Perry (Ref. 6] develops his Lattice modeling algorithms. 
Equation (2.30) can also be recognized as the solution to 


the Yule-Walker or Normal equations [Ref. 26]. 


2. Moving Average Form (D. = Q) 


Referring once again to equation (2.7) and setting 


m@emeoefficients b. = 0 for all i we have case (2), 


Vase ee) (ei) 6+ v(k) (251) 
i=0 
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Since v(k) is by definition noise associated with the 
measurement, we can combine y(k) and v(k) such that, 


z(k) = y(k) - v(k) = a,u(k-i). (2-62) 


aS 
© 


Equation (2.32) simply states that the present measurement 
is a linear combination of past and present inputs or by 
definition, a Moving Average (MA) process [Ref. 25,33]. The 
Adaptive Kalman Identifier estimate for the coefficients of 


the MA process is as before, 


a(k+1|k) = {FF =K(K)H(k)] a(k+1|k) + K(k) 2(k) (238 Srey) 
ut aE -1 

K(k) = P(k|k-1) H’(k) [H(k)P(k|k-1) H™(k) + Rl (272 315) 

P(k+1|k) = P(k|k-1) -K(k)H(k)P(k|k-1) +Q. (252) 


Using the alternate forms of the above equations we arrive 


ac, 
a(k+1|k) = [P(k+1|k) (P(k |k-1) ] 7] a(k |k-1) 
T a 
+ (P(k+1|k) H?(k) R71] z(k) (2.34a) 
Sik) = ({P(k|k-1)] + +H %(k) + Htk)]?. (2.34b) 


Arguing as we did for the autoregressive case, no a-priorl 
knowledge about the initial values of the moving average 


coefficients, implies that, 
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reo iles = 0. (2.35) 


And even if after we processed one measurement we still 
admitted no knowledge as to the accuracy of the previous 


estimate, we imply that, 


4 (0s ei) |) gee ae Oe (2.36) 


Substituting these implications into equations (2.34) our 


estimate becomes, 


i 


a(k+l|k) = ([H(k) R72 a7 (k)]7+ 


H(k) R°- 2(k) (2.37) 


the by now familiar weighted least squares estimate [Ref. 28, 
29,30). If once again we allow al to be unity, let k tend 

toward infinity and take the expectation of equation (2.37), 
one arrives at the discrete form of the Wiener-Hopf equation 


[Ref. 11], namely 


“ 2 -1 
a(k+l|k) = RU (k) x(k) (2.38) 


where it can easily be shown that, 


Ray 6?) Raut) Ray om) 
Sep Oa Ste we Ray (tom) 
Ray *) = ; ; : C23 9) 
Roy (7m) eee Ry ' 9? 
and, 


Zo 





uy 
gg 

‘ea = : : (220) 
Nyep 


Meeaclton (2.38) provided another point of departure from 
which Perry [Ref. 6] develops the MA Lattice modeling algorithms. 
Appendix A develops the Wiener-Hopf equation otherwise known 
as the all-zero model from yet another approach. 

The all-zero model is fundamental to Widrow least 
mean square (LMS) adaptive filters [Ref. 11] and linear pre- 
diction theory [Ref. 16]. 

pe Autoregressive-Moving Average Form 

Returning to the alternate form (equations 2.20) of 
the Kalman Filter equations (2.19) the cevelopment of the 
Autoregressive Moving Average (ARMA) Adaptive Kalman Iden- 


tifier follows. The estimate for the as and b. coefficients 


1s, 
a (k+1|k) a(k|k-1) 
a. SGP (k (kel) ) 1], |-=------- 
b(k+1|k) b(k|k-1) 
T ail 
+ [P(k+l|k) H™(k) R ~“]y(k) CoE 
Beeiik) = [(P(k{k-1)) 1 +H2 (x) Ro? H(k)17+ . (2.41b) 
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Progressing in the same manner as in the previous two cases, 


we assume, 


(= (0 |) een (240 )a)) 
-l1 ._ 
PeCelk—19a) y= 0 (2.42b) 
The estimate then becomes, 
a (k+1[k) 
a Seer H(i eH (hk) Roy Ck) - (2.43) 
b(k+1 |k) 
Taking equation (2.43) a step further by letting aoe be 


unity, letting k approach infinity and taking the expectation 


we have, 
A -l 
a(k+1/k) Ray ‘*) Beagy USE) Seg ee 
= - woe ------ (2.44) 
a i 
e(k+1)k -R k-1 R k R k 
D |k) uy ) yy ' ) vy ‘ ) 
Note that the time varying measurement vector, H(k), is of 
the form 
H(k) = (u(k) u(k-1) ... u(n-m) | -y(k-1) -y(k-2) ... -y(k-n)] 
(2.45) 


Shown in detail, equation (2.44) has the characteristic form 
(equation 2.46, next page). It is the Toeplitz and symmetric 
Nature of equation (2.46) that is exploited by the Levinson 


[Ref. 34] and Lattice [Ref. 35] algorithms. 
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cr Sn nt ae ee ee, f- ee ee ne ee ee = | aves 
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The Adaptive Kalman ARMA modeling technique can be 
eae Visualized by a block diagram. Referring to Figure 2.1, 
the unknown system is excited by a white, zero mean, gaussian 
noise sequence of sample values from a random process. The 
input and its associated output are then passed through M 
and N delays respectively in serial form. The parallel inputs 
(M+l1l) and the outputs (N) are concatenated to form the 
measurement vector H(k) which is represented by equation 
(2.13). The inputs N and M are selected a-priori as the 
model orders for the Autoregressive and Moving Average pro- 
cesses one desires to identify. It can be easily seen that 
the remainder of the figure simply implements equation (2.15a). 
The outputs of the Adaptive Kalman Identifier are an estimate 


of the coefficients, 


and a one step prediction, y(kj|k). 


©. OBSERVATIONS 

In this section it has been shown that a direct connection 
can be established between the Adaptive Kalman Identifier and 
the Yule-Walker equations associated with an AR process. 
Secondly, the steady state Adaptive Kalman Identifier closely 
resembles the discrete Weiner-Hopf equation associated with 
the MA process when the measurement matrix, H(k), is of the 


BOM , 
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Beene kK—-1) u(k-2) ... u(k=m) J 


And, thirdly under the same assumptions as in the previous 
two cases, the Adaptive Kalman ARMA Identifier is similar to 


the form Perry [Ref. 6] exploits using Lattice modeling. 
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IIIT. COMPARISON BETWEEN THE ADAPTIVE MA KALMAN 
TDENTIFIER AND THE WIDROW LMS ADAPTIVE FILTER 

A. PRELIMINARIES 

It is instructive to investigate the similarities between 
the Adaptive Kalman Identifier and the Widrow LMS Adaptive 
Filter when the concepts are applied to system identifica- 
tion. Basically, the LMS algorithm implementation of system 
identification considers a block diagram as is shown in 
Mere 3.1. The output, y(k), of the Adaptive filter 1s simply 
a weighted linear combination of the past and present known 
inputs. The same input is fed into both the unknown system 
to be identified and the adaptive filter. The output of the 
unknown system is designated the desired response, d(k), from 
which an error signal, e(k), is derived. The error signal, 
e(k), provides the criterion through which the weights, Wey 
are adjusted such that the error, ¢(k+tl), is driven toward 
its minimum. A more detailed analysis of the operation of 
the LMS algorithm can be found in Ref. ll. 

The weights w, are aewusted from eime Step to Eine 


step in the following manner, 
w(k+1l) = w(k) - 2k. e(k) X(k) C3215 


where we define the following quantities, 
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Wo (k) 
w, (Kk) 
Oo : (222) 
Wik) 
the weight vector at time step k, 
Gk) 
meek =>) 
Atk) = , ( 35s) 
u (k-m) 


the input signal vector at time step k, 


ea kK) 


Gites] vik) (3.4) 
y(k) = w'(k) X(k) els 


mmemerror, €(k), and the output, y(k), at time step k, and, 
Kos a scalar constant controlling the rate of convergence and 


the stability of the adaptive filter. 


Be THE COMPARISON 


Meester ctucing (3.4) and (3.5) into (3.1) gives, 


w(kK+1) = w(k) -2k_ X(k) [d(k) -w' (k)X(k)]. (3.6) 
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Recalling the MA form of the Adaptive Kalman Identifier, 


equation (2.33a), rewritten here for convenience as 
a(k+l|k) = ead leer Ws) iad) - (Kk) a(k|k-1) 1, a7) 


one can make the following associations: 


a(k+1|k) (= w(k41) (3. 8a) 
a(k|k-1) (= w(k) (3. 8b) 
z (kK) (= a(k) (3. 8c) 
T 

ne (k) ex (x) i2eed) 

K (k) = -2k, X(j) (3.8e) 


Recall that for the MA form of the AKI, the measurement 
vector H(k) represents a vector of past and present inputs. 


Namely, 
Teo uik) u(k=1) ... u(k-m)]j . ea 


Therefore the associations (3.8a)-(3.8d) are straightforward. 
However, it iS not so clear as to what is meant by the 
association (3.8e). Digressing a moment to present an equi- 


valent expression [Ref. 27] for the Kalman gain, K(k), 


K(k) = P(klk) Ho (k) rR? (3.10) 
and substituting (3.10) into the association (3.8e) and enter- 


taining the conjecture that the quantities are equivalent 


a9 





under certain conditions, we have: 


P(k|k) H’(k) Ro = -2k, X(k) . (3.11) 
The conditions alluded to are (1) that the Adaptive Kalman 
Identifier is in steady state and (2) that the statistics of 
the input forcing function are stationary. Denoting the steady 
state error covariance, Ba ik) ase .; equatiena(3.11) can 


be solved for Koy 


KI = - = Eee: (3222) 


Invoking the entire Kalman gain history in its more popular 
form of equation (2.33b) and equating it to the Widrow gain 


(3.8e), we instead arrive at, 


-2k_X(k) = P(k[k-1) H7(k) [H(k) P(k|k-1) H™(k) +R] 
(3.13.4) 
T 
Hl(k) = X(k) . (3.13b) 
Solving for the convergence factor, Kos We wobtcain, 
“2k X(k)X"(k) = P(k[k-1) X(k) [X"(k)P(k|K-1)X(k) +R] 7X" (k) 


cg Gel 


-2k_I = P(k|k-1)X(k) (x7 (k)P(k|k-1) X(k) +RI 7X7 (k) (x(k) x7 (k) J 7+. 


But [x (k) P(k|k-1)X(k) +R]? So a siecle... ENererore , 


P (k|k-1) 
klo= (3.14) 
7 X°(k) P(k|k-1) X(k) +R 
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Griffiths [Ref. 22] suggests that the convergence factor Kes 
which he denotes as u(n), should be chosen such that, 


u(n) = ———— (32115) 
L Oo, (n) 


where 0 < a < 1 is a normalized adaptive stepsize parameter 
and the term 5° (n) is an estimate of the input power level 
which may be computed using a geometrically-weighted average 


for an L weight adaptive filter, 


A2 i oc oo OL ee 
o.(n) = (1 ) oO. (n 1) + r <o-(0)) : (36) 


Comparing equation (3.14) and (3.15), one observes that if 

P(k k-1) is equal to the identity matrix then we obtain a 

term proportional to the input power. Further if measurement 
noise is considered negligible, R ~- 0.0, then 3.14 and 3.15 

are essentially equivalent. The salient point to make, however, 


is that the choice of the convergence constant, k in the 


Ss! 
LMS algorithm provides no clue as to how to deal with measure- 
ment noise whereas the Adaptive MA Kalman Identifier explicitly 


incorporates measurement error into its algorithm. 


me OboERVATIONS 

It has been known that the LMS algorithm was suboptimal 
Since the actual gradient of equation (3.1) is replaced by 
the approximation, 2€(k) x(k), [Ref. 24]. However, the degree 
memouooptimality is difficult to quantify. In previous 


linear prediction research no mention was made regarding the 
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role of measurement noise. Equation (3.14) gives us some 
insight into wiser selections of the rate of convergence con- 
stant, Ker Since it takes into account the effects of measure- 
ment noise. Further, it appears that the LMS adaptive filter 


is a special case of the Adaptive MA Kalman Identifier. 
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ITV. COMPARISON BETWEEN THE ADAPTIVE ARMA KALMAN IDENTIFIER 
ANDetaE ADAPTIVE RECURSIVE LMS FILTER 

This chapter endeavors to investigate the similarities and 
differences between the Adaptive ARMA Kalman Identifier and 
the Adaptive Recursive LMS Filter [Ref. 18]. LMS filters 
[Ref. ll] have enjoyed much popularity in the recent past due 
to their ease of implementation, simple unimodal algorithm, 
robustness and ability to "adapt" to the unknown statistics 
of the signal environment. Being transversal in nature, they 
have a finite impulse response being able to produce only 
zeros in the input/output transfer function. One may, how- 


ever, decide to model the transfer function, 


i 


H, (2) = 5 ena (A sale} 
- .92 


using a transversal filter, only to realize that a large 
number of delays are required in order to arrive at an 
adequate approximation. The germane point which is being 

made is that one weighted feedback tap can realize an infinite 
string of feed-forward coefficients. Moreover, it 1s very 
desirous to adjust the feedback weights adaptively in some 


optimal fashion to the statistics of the signal environment. 


A. PRELIMINARIES 

The approach taken by Feintuch [{Ref. 18] is patterned 
after the analysis first presented by Widrow [Ref. ll]. A 
Summary of Feintuch's derivation will be presented here with 


emphasis on its application to ARMA modeling. 
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Recall the general ARMA equation (2.1b), 


= 


n 
ee) = ) a.u({k-i) + } b,y (k~i) (4.2a) 
i=o ?+ Teoh 


which can be rewritten in vector notation, 


y{k) = au x by (4.2b) 
where, 
a" = [ay ay a! (asa) 
oe Saya b ] (4.35) 
— rapes n ? 
T 
Wao yk lL) yYtK=2) =... y(k=n) } (4236) 
T 
ee i (cet k=l) .. . u(k=m) | (453) 


Given that (4.2) is the assumed mathematical description of 
PeemunkKnown system where u and y, the input and output data 
Meeeemces, are known, then by solving (4.2) for a and b, 
identification of the unknown system can be made. The LMS 
algorithms, in general, employ a "desired" signal d(k) with 
which to "train" the adaptive filter. If the desired signal 
waS assumed to be the response of some unknown system to a 
known input signal, then the algorithm presented by Feintuch 
[Ref. 18] can be used as a means by which to identify the 
system parameters. 

The first step in the Recursive LMS derivation is to form 


the error between the desired signal, d(k), and the output 
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of the filter, y(k), 


e(k) = d(k) - y(k) = 4(k) = alu - b’y (4.4) 


Forming the square of (4.4) and taking the expectation, 


we have the mean square error representation of the filter, 


Bfe“(k)} = E{d*(k)} +a" R a+b Be: <>. 
5 RS Ase Ray + 2a* RayP (4.5) 
where, 

R (kK) = Etu(k) u’(k)} (4. 6a) 
Ryy(K-l) = Ety(k-1) yi (k-1)} (4.6b) 
Ra, (k) = E{d(k) u(k)} (4.6c) 
Ray(k) = E(a(k) y(x)) (4.64) 
Ryy(k-1) = E{u(k-1) y(k-1)} . (4.6e) 


It 1S assumed that the statistics (4.6) are constants allowing 
the gradient of (4.5) to be taken with respect to the a; and 
D.. This assumption does not stretch the theory since in 
practice one uses an input test signal with stationary char- 
acteristics and if the unknown system's output process is 

not stationary, then no identification of the system parameters 


Can be made. The respective gradients are, 
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3 2 = fa 
x, Ete (k)}] = 2R oy (Ke ge Sa (ara) 


9 2 2 ab 
ap Ete (k) 3] aR. (k) b 2 y (Kk) +2R uy (K) a (Qapey, 


Ra 


Since the second order statistics are assumed to be known, 


Semetions (4.7) can be solved for a and b by Setting the 


Gradients equal to zero. In matrix partitioned form we have, 
k k 
Eu ) ! uy ) a Ray (*) 
a ———— | i cee 225 = 2 : (4.8) 
R?  (k) | (k)| | b R,_ (k) 
uy aay, = dy 


At this point we digress momentarily to mention that Johnson 
and Larimore [Ref. 19] and Widrow and McCool [Ref. 20] agree 
with Feintuch's derivation. The following steps, however, 
are controversial [Ref. 19,20]. 

Feintuch continues and states that in general, the statis- 
tics involved in (4.8) are not available a-priori. One method 
of estimating the statistics is to make the filter adaptive 
in an iterative fashion using the method of steepest descent. 


The method of steepest descent employs an algorithm of the 


ZOrm , 
Mey = alk) + k. *-[Rle*(k)}] Ceca 
— — ada . 
_ P) 2 
b(k+1) = b(k) + k, SplB(e“(k)}1. (4. 9b) 


The gradient involved in (4.9) can be approximated using the 


techniques outlined in [Ref. 11]. The iterative LMS algorithm 
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Meeeectinating the coefficients a and b is, 


a(k+1) a(k) - 2k e(k) u(k) (4.10a) 


b (k+1) b(k) - 2k,e(k) y(k). (4.10b) 


Figure 4.1 is a block diagram implementation of the Adaptive 
Recursive LMS Filter algorithm as applied to system identifi- 
cation. The input signal, u(k), to both the unknown system 
and the adaptive filter is a zero mean, white gaussian noise 
sequence of samples from a random process. The output of 

the unknown system is d(k) which is used as the training 
eeonal for the filter. The output of the adaptive filter, 
y(k), 1s compared with the desired signal, d(k), from which 
an error Signal is derived. The error signal, e(k), is then 
used in the LMS algorithm, equation (4.10), to iteratively 
Meee tne a and b coefficients. Theoretically, when the 
coefficients of the zeros-producing adaptive filter and the 
poles-producing adaptive filter correspond exactly with those 
of the unknown system, the error sequence should be zero and 
the unknown system is modeled by the Adaptive Recursive LMS 


Filter. The convergence constants, Ko and k are chosen 


Bb” 
Meee rt as are the initial values, a(0) and b(0), for the 


coefficients. 


B. THE COMPARISON 
One can begin the comparison between the Adaptive Kalman 


Identifier and the Adaptive Recursive LMS filter by noting 
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the similarities between the equations (2.43) and (4.8), 


repeated here for convenience, 


| 
a(k+1|k) Ray (*) | ee) Ruy (Kk) 
— = SS a a ————— (43a) 
k k - k-1 R k k 
b(k+1/|k) R ae “| ) wan ) 
A, = B, Cy (2.43b) 
| 
Ray bk) | Ruy (koL) a Ray (*) 
a : a = = ------ (4.8a) 
T 
k-1)! R_ (k Bua bs 
R ny | ue ) b au) ) 
3 R_ (k) | 2 (Oren) | ve Tne 
a uu | uy du 
= a : m--------/|} #| ------ (4.8b) 
A T 
b R k-1) |! R_  (k R. (k 
b ms - ar ) aug ) 
A, = B. C., (4.8c) 


ime upper left partitions of By and B. are obviously identi- 
cal, in that the AKI and the Adaptive Recursive LMS both 
implement instantaneous estimates of the input covariance 
Eunections, Roy bk) - Tie ewer right partitions of By and B., 

are Similar with one subtle difference. Both are instantaneous 
covariance functions; however, By employs the covariance, 
Bey) of the actual data whereas B. uses, in fact, previ- 

Ous filter estimates to compute the output covariance, aaa 
more properly denoted by Se Sd The negative sign of the 
Upper right and lower left partitions in By, in comparison 

with B is a result of the initial definition for the recursive 


pi’ 
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weights, b,. That Pp lemerains ser nrumct1on fOr the general 


ARMA equation can be written as, 


th -i 
Be Paez 
Hapmal2) = —— 45 (4.9a) 
1- ) b Zee 
i=1 7? 
or, 
=). 
ag + Pe. a. Z 
H" . ema ‘2? = = 7 (4°95) 
Te Sayers 
; o% 
ese 


The second and more important observation is that cross- 
correlations employed by equation (4.8b) use the instantane- 
@uemvalues for the past estimates of the outputs while (2.43a) 
uses the actual past values for y(k). Similarly, comparing 

Cy with C., SMeweinass thac tne lower partition of C., uses the 
cross-correlation between the desired signal and the past 
estimates of the adaptive recursive LMS filter. 


Equation (4.8a) can convey the above information more 


Clearly if it is instead written as, 


x } 
j 
a aul) TE Rok) || Rg (0 
— a op = 8&3 a == aE aD amp ape ape 1 ap epempemepemw ep emmamm —§|— £iéi$|$ of om a 2 a= a= (4210) 
- Al 
b R A k= j AA k A k 


The convergence factors, ie and k can be compared to 


Jove 
the steady state Adaptive Kalman Identifier gains by making 
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the identical analogies (3.8) as was done in the previous 
chapter. One also observes that since the output data that 

is processed by the Adaptive Recursive LMS filter are in 
actuality filter estimates, a more correct version of Feintuch's 


Seroginal algorithm (4.10) is written as, 


a(k+1) eu) i 2k. eves) uk) (4.lla) 


b(k) - 2k, e(k) y(k) . (4.11b) 


b(k+1) : 


In agumented form, equation (4.11) is, 


a(k+1) a(k) 2k. u(k) 

------ = | ----]| - | -------- |[a(k) -a(k) u(k) - b(k) y(k)] 

b(k+1) b (k) 2k, y (k) 

a(k+1) a (k) 2k_ U(k) a (k) 

----- = | ----| - | --------]}]a(k) - [u(k) y(k)] |---- 

b(k+1) b (k) 2k, y (x) b (k) 
(4.12) 


Recalling equation (2.15a) in a slightly different form we 


have, 
a(k+1|k) a(k{k-1) a(k|k-1) 
_—------ = |--------| + K(k) |z(k) -H(x) |-------- 
b(k+1 |k) b(k |k-1) b(k|k-1) 


CaaS ) 


It is a simple matter to explore the similarities between 


(4.12) and (4.13) and make the following associations, 


ie 





AKI Adaptive Recursive LMS 


a(k+1|k) a(k+1) 
_—-—..- eS ————— (22a) 
b(k+1|k) bck) 
a (k|k-1) a (k) 
a __ (E59 peas (4.14b) 
b (k|k-1) b (k) 
z (k) uae) (4.14c) 
meee) 6 = =SCs [u(k) ... =) [u(k) y(k)] (4.14d) 
Goat) y (ka1)} 
Vike) 
Pras) 
K (Kk) C=) | eee (4.14e) 
2k Yy(k) 


The associations (4.14d) and (4.14e) are the major differences 
between the AKI and the Adaptive Recursive LMS. It has been 
Shown that the mean square error surface employing not only 
estimated feed forward coefficients but also estimated feed- 
back (recursive) coefficients is in general multimodal [Ref. 
36]. Hence, the Adaptive Recursive LMS algorithm does not 
minimize the mean square error [Ref. 19,20], and in general 
the gradient algorithm in this case does not seek the global 


minimum. 
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eee OBSERVATIONS 

The complex nature of estimating recursive coefficients 
memyet a formidable problem [Ref. 19,20]. It seems that 
Feintuch's algorithm is successful due to the a-priori knowledge 
of the minimal order generating the desired signal. This 
knowledge is used in setting the order of the adaptive 


Mecursive filter. 


ae 





V. SYNTHETIC DATA GENERATING PLANTS 


The Adaptive ARMA Kalman Identifier was tested using 
computer derived data from several models. The input/output 
data was collected by driving a known plant with a zero mean, 
unit variance, white gaussian noise sequence of samples. 

The AKI algorithm was repeatedly tested using data from 


progressively more complex models. 


A. AUTOREGRESSIVE-MOVING AVERAGE DATA 
Autoregressive-Moving Average (ARMA) data was generated 


using an equation of the form, 


y(k) = byy(k-1) + boy (k-2) gfe eaeetae ects bY (k-n) 


+ a.zu(k) + ... + aja (k-m) (Seu) 


0 
Taking the Z-transform of (5.1) we have, 


2) = byz ¥ (z) ~ poz °¥ (z) + ... + b 2 ¥(z) 


eevee 2. + a2 U(2) (Se) 


where Y(z)(U(z)) represents the Z-transform of the output 
-n -n 
meieut) and z Ov (2) (2 


°U(2z)) represents the Z-transform 
of the output (invut) delayed No time steps. Equation (5.1) 
has already been defined as the general form of the ARMA 


process provided that the sequence u(k-i) i = 0,1,... comes 


from a gaussian random process [Ref. 25]. The ratio, 
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-l 
a. +ta.,Z > ee = ee 4 
Ee a ike a (5.3) 
ri Z = ee © ee 

1 n 





can then be readily recognized as the general transfer func- 
Seem (Ref. 37] of a digital plant. That is, the u(k) is the 
observed input and the y(kK) is the perfectly measured output. 
Consider that the measurements of the output include some ran- 
dom error, v(k), due to the inaccuracy of the meaSuring instru- 
ments, then the resulting situation is as depicted in Figure 
5.1. A practical, judicious and reasonable description of 
the measurement error is that this be a stochastic process 
whose distribution is zero mean gaussian. The noisy measure- 
ment, z(K), 1s then the output, y(k), plus the measurement 
mmmor, v(k). 

There are three cases of interest: Case l: b. = 0 for 
Seeeecase 2; a. = 0, 1 = 1,2,..-.,m; Case 3: a. #0, 


1 1 
b. mo, cor all 1. 


H Moving Average Data (Case 1) 
Case one 1S readily recognized as the all zero plant 
which produces moving average data. A simple second order 


plant was defined where a, = 1.0, a, = 2.0 and a. = 3.0. The 


0 IL Z 
gaussian noise sequences, u(k) and v(k), were obtained using 
the general purpose IMSL subroutine, GGNML, provided by the 
computer center at the Naval Postgraduate School. 

2. Autoregressive Data (Case 2) 


Case two results in the all pole plant producing auto- 


regressive data. The coefficients were selected such that the 


De 





meememot Figure 5.1 produced stable data. In precise terms, 


the poles of the function, 


= Jk 
o- 4,26 ,-1,9 2,1 ,-3 poe) 


DE waa OA 


1 


were located completely within the unit circle in the Z-plane 
ensuring a stable plant [Ref. 37,38]. The noise sequences 


were obtained as before. 


v (k) 


u(k) 





Figure 5.1. ARMA Digital Plant 


3. Autoregressive-Moving Average Data (Case 3) 

Case three follows from a logical combination of the 
two previous cases. That is, the transfer function, H(z), now 
has both zeros and poles. In order to compare with previous 
work done in the area of system identification [Ref. 6], one 
Of Perry's models [Ref. 6] was used. Namely, the transfer 


function H(z) used for this case was, 


oul =e 
2) «(= IES AUS eA eee sly. z (3.5) 


ene erections’ = 694902 °° + .407452 
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The input, u(k), was obtained as in the previous two cases 
using the IMSL subroutine GGNML set to unit variance and zero 
mean. The output measurement noise sequence, however, was 
set to a variance of .0001. The basic software, program ARMA, 
which generated the different data is included as Appendix B. 
Basically, the coefficients A1l-Al1l0 change the character of the 
general ARMA equation (5.1) which can be selected to produce 
the desired plant data. The input/output data is then written 
onto a disk file for subsequent analysis. 

In order to compare some of the findings in this the- 
Sis with the results previously obtained by Feintuch [Ref. 18], 


Mme transfer function, 


-1 
fie) = 0.05 -0.40z ae 


iL ol Sane oO sen 


was also used as a source of synthetic data. Feintuch used 
(5.6) as a source of data with which to analyze the operation 


of his Adaptive Recursive LMS Filter. 


Bee PHASE LOCKED LOOP DATA 

The second class of data used to exercise the Adaptive 
ARMA Kalman Identifier was derived from a computer simulated 
phase locked loop (PLL) developed by Romeo [Ref. 7]. The 
basic PLL algorithm, however, implemented a forward Euler 
integration scheme. 

Briefly, in block diagram form, the phase tracking 
Characteristics of the PLL in the frequency domain can be 


depicted as in Figure 5.2. 


a 








Figure 5.2. PLL Block Diagram 


Where: 


A is the input voltage 
K is the loop gain 


Fis) is the filter transfer function 


For most applications, it is assumed that o(t) << 7/2 allowing 
meento Make the approximation, Sina *a. This is the 

linear mode of operation. Romeo [Ref. 7] chose the filter 
Characteristics F(s) = 1+ K/S and the same characteristics 
were used in this thesis in order to compare results. The 
Parameters of the overall system were adjusted to obtain 
mecempang coefficient, ¢, of 0.3 and a natural frequency, Wie 
of 3.33 rad/sec. This resulted in a step response overshoot 
memcmout 46.7% at t ~ .75 sec. Solving for the time domain 


step response analytically in the linear region one obtains, 


ae 


mee) 6 = 1 + 1.048e ~ sin(3.178t -1.266) (5.7) 


The final block diagram of the PLL system is shown in Figure 5.3. 
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Y(s) 





Beoaure 5.3. Final PLL Block Diagram 


Normalized step responses for several step magnitudes 
uSing Digital Simulation Language (DSL) were obtained. Figure 
5.4 shows the step responses of several step inputs in the 
meme 0 < u(t) < 30. It is obvious that the PLL is in its 
approximate linear region and does not exhibit any discernible 
Spoetortion. Figure 5.5 is a repeat of the above test; how- 
Mueeecne range of step inputs are 30 < u(t) < 170. It is 
apparent from Figure 5.5 that the sin (-:) nonlinearity begins 
to have some effect on the operation of the PLL in that the 
amplitude of the responses are reduced and delayed. 

The DSL PLL system, however, was not used as a source of 
Synthetic data because of the problem of nonstationary statis- 
tics of the input signal using the random sequence generators 
aVaillable under DSL already discussed in Ref. 7. These 
Simulations were nevertheless used as a basis for comparison 
of the operation of the discrete PLL developed in Reference 7 
but modified to implement a forward Euler integration scheme 


rei 1S used in this thesis. 


a 





MmmegoAROER PLL -- FRAPEZOIOAL INTEGRATICGN 
MyemMALIZED STEP RESPONSE COO / 1/4 
00 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 

Poerec= 1.00 UNI TS/INCH RUN NQ@. i 

[eenee= 0.20 UNITS/ INCH i Ieee) 


Figure 5.4. 


PLL Step Response, Linear Region 


60 


C7 207 30 








SN GRADER PLL -- TRAPEZOIUQURL INTEGRATION 
Meili ZE® SIEP RESPONSE 1730/60/90/120/150/1 7/70 
“ | 

i ¥ MW S 

\W 
00 1.00 2.00 3.00 4.00 5.00 §.00 7.00 8.00 9.00 
Mmeenee= 1.00 UNITS/ INCH AUN: NG 1 
MmeenLE= 0.20 UNITS/INCH mae] Nias | 


Figure 5.5. PLL Step Response, Non-Linear Region 
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The forward Euler integrator takes the form, 


-l 
toa 2 _ a tol (5.8) 
Ss -l 
1 -2z 


The block diagram of Figure 5.3 changes character to appear 


Memin Figure 5.6. 





Figure 5.6. Discrete PLL Block Diagram 


mepeying Mason's gain rule to the block diagram of Figure 5.6 


in the linear region we arrive at, 





-l -2 
. een O02 027 - ee (5.9) 
dh IEP f6' (0102 eee Saul OZ 
the linearized PLL transfer function. The discrete PLL sys- 


tem has zeros at (0.0,.9445) and a complex conjugate pole 
Metre at (0.990 + 3.03178). 

The discrete PLL system was implemented uSing a Fortran 
program on the IBM 370 in double precision. The source code 


for the discrete PLL is included as Appendix C. The discrete 
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PLL was tested uSing various magnitude step inputs. Results 
were Similar to those from the DSL simulation. That is, the 
Me crete PLL exhibited the same approximate overshoot, rise 
time, natural frequency and nearly identical time of max over- 
shoot in the linear region. To implement the PLL using forward 
Euler integration, the block diagram of Figure 5.6 was recon- 
figured to appear as Figure 5.7. Note that Figure 5.7 does 


not employ any delay free loops, hence it 1s easily programmable. 


= 
1-.9445z y (k) 





Figure 5.7. Programmable PLL Implementing 
Forward Euler Integration 


Figure 5.8 summarizes the results of the discrete PLL simula- 
ten for normalized step inputs in the range 0 < u(t) < 30. 
Similarly, the discrete PLL was tested in the nonlinear region 
meeagenormalized step inputs in the range 30 < u(t) < 170. 

The discrete PLL displayed the same characteristics of de- 
creased output amplitude and delay which was first observed 

in the DSL simulation. Figure 5.9 summarizes the test results 


of the discrete PLL in the nonlinear region. 
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Figure 5.8. Discrete PLL Step Response, Linear Region 
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Discrete PLL Step Response, 
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It was felt that the PLL was adequately modeled for use 
as a source of data to be subsequently analyzed by the Adap- 


tive Kalman Identifier. 
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VI. MODEL IDENTIFICATION SOFTWARE 


Three basic Fortran IV programs for use on the IBM 370 

were developed: 

(1) ADAPTSN: Adaptive Kalman Identifier 

(2) LMS: Adaptive Least Mean Square Filter 

(3) LMSR: Adaptive RecurSive Least Mean Square Filter 
All three programs used double precision arithmetic to minimize 
mm@enetrects of truncation error, limit cycles and roundoff 
error. In all cases the orders of the autoregressive and 
moving average processes were read in unformatted form from 
a disk file along with the input/output data of the unknown 
system to be analyzed. An attempt was made to minimize memory 


storage, but not at the expense of program flow and clarity. 


A. ADAPTIVE KALMAN IDENTIFIER 

Program ADAPTSN is a versatile Fortran IV software program 
which implements equations (2.15) to identify the coefficients 
of an MA, AR, or ARMA process. Additionally, at each iteration 
the poles and zeros of the evolving transfer function are com- 
puted. ADAPTSN can also be used to perform regression analysis 
On non-linear terms as discussed in Chapter VII. The versa- 
tility of the AKI lies in the various options which are avail- 
able to the user. By properly selecting the flags N, M, and 
NL, the user can perform regression analysis on data whose 
combined order (N+M+l) is less than or equal to 20. Table 6.1 


lists the ascribed meanings of the various flags. 
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mame 6. Il 


AKI FLAG MEANINGS 


N: the order of the autoregressive process 


M: one greater than the order of the moving average 
process 


NL = 0: AKI is used to identify the coefficients associated 
with the general ARMA equation (2.1b) (linear 
regression) 


lI 
ra 


NL AKI is used to identify the coefficients of the 
linear general ARMA equation and the weighting 


coefficient associated with wu” (k-i) 


NL = 2: AKI is used as in NL = 0,1 and additionally iden- 
tifies the weighting coefficient associated with 
y3(k-i) 

NL = 3: AKI is used as in NL = 0,1,2 and additionally iden- 

tifies the weighting coefficient associated with 

u2 (k-i)y (k-i) 

AKI is used as in NL = 0,1,2,3, and additionally 

identifies the weighting coefficient associated with 

y2(i-i)u(k-i) 


iI 
> 
ee 


NL 


NL = 5: AKI is implemented to analyze time series and 
identify the ARMA coefficients associated with 
the Box-Jenkins model [Ref. 26]. 


mie Option given by NL = 5 was not extensively tested and as 
such only limited results are available. Table 6.2 outlines 
the allowable flag combinations which will produce a valid 
analysis of the data. 

For purposes of this thesis, the non-linear (NL) options 
were configured to implement the expansion terms associated 
With the Taylor series expansion of the sine function. This, 
however, can be changed at the user's discretion to implement 
any other series expansion by inserting the appropriate Fortran 


Statements at the proper location in the AKI program. 
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Table 6.2 


VALID FLAG COMBINATIONS 


OPTION (NL) MAX COMBINED ORDER 
Q Nee ie 2 0 
ak N+ 2M < 20; N,M # 0 
2 2N + 2M < 20; N,M # 0 
3 3N + 2M < 20; N,M # 0 
4 4N + 2M < 20; N,M # 0 
5 Nee Ms 20 


The overall structure of the AKI program uses subroutine 
calls to compute the various quantities necessary for the 
eventual computation of the system ARMA coefficients. These 
Subroutines are: GAIN, RECKON, PRINT, NEXT. A brief des- 
Ccription of the function each performs is given at the 
beginning of each subroutine. Figure 6.1 describes the 
general program flow of the AKI. The source code for program 


ADAPTSN is included as Appendix D. 


ee ADAPTIVE LEAST MEAN SQUARE FILTER 

Program LMS uses the equations developed in Chapter III 
to compute the coefficients (or weights) associated with a 
moving average process. The LMS filter program is capable of 
computing up to 12 MA coefficients; or equivalently, is capable 
of identifying an eleventh order moving average process. The 
general program flow of the LMS filter is presented in 
Figure 6.2 and the annotated source code is included as Appen- 


max EF. 
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Figure 6.1. AKI General Program Flow 
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Figure 6.2. LMS Filter Program Flow 





meen DAPTIVE RECURSIVE LEAST MEAN SQUARE FILTER 

Program LMSR realizes the algorithm proposed by Feintuch 
and recapped in Chapter IV. It was designed to handle a 
combined AR and MA order of ll. That is, the program can be 

y 

used to estimate m+ 1= M coefficients associated with the 
MA process and n = N coefficients associated with the AR 
Meeeess SUCH that N+M< 12. The general program flow for 
the Adaptive Recursive LMS filter 1s shown in Figure 6.3 and 


the source code is included as Appendix F. 
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Mit. NON-=ETNEAR IDENTIFICATION 


The Adaptive Kalman Identifier can be modified to estimate 
the weighting coefficients associated with higher order terms 
produced by some non-linearity within the unknown system. In 
general, much has to be known about the non-linearity such 
that the functional description chosen for it is a close 
approximation to the effects it causes. In this thesis one 
possible application of the Adaptive Kalman Identifier toward 
identifying a system with a known non-linearity is explored. 
No attempt has been made to present an extensive treatment of 
non-linear analysis techniques. 

The approach taken has been previously explored by Parker 
(Ref. 8]. A special case of the generalized non-linear ARMA 


model [Ref. 8], 


y(k) =_ J SU TL OMS ea ) | ) a(i,,i,)u(k-i,)u(k-i,) +... 
1,=0 We=0 1 2=0 
i iL 2 
-_ lL ; patty. si) Goa, )) 2 Sega 
= m 
+  ) bUj,y(K-5,) + YY BG yr dg)y (k-G,)y¥(k-55) +... 
a 3,71 j5=1 
2, ee pega, )yte-j)ytk-3,) ..-¥(k-3,) 
a0 2 
. : L | D c(i, ,j,)ulk-i,)y(k-j,) + - a ; 7 .; 
4 afl iL m Jy 
| ) c(i, +---d,5,---5,)ulk-i,)-..ulk-i,) y(k-j,)..-y(k-j,) 
Drs 
Cay) 


es 





is used to model the input/output relationships, non-linear 
transfer function, of the phase locked loop. The non-linear 
element, the sine function, of the PLL can be replaced by 


a Taylor power series expansion, 


pen(x) = xX - + x? ~ = x> a te (ce ) 
Other expansions can be used, e.g., Legendre polynomials, 
Volterra series, ..., etc.; however, only the Taylor expan- 
Sion was investigated in this thesis. Practical implementa- 


meen Of {/.2) suggests truncation at the third order term. 


The sine function is therefore approximated as 


3 


Sin(x) = x - = xX 7 ea) 


Substituting (7.3) into the block diagram of Figure (5.6) we 


have Figure 7.1. 






I 





-k 
BZ u(z) + ae 


Figure 7.1. PLL, Third Order Taylor Series 
Approximation 
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The discrete time domain equations at the different nodes are: 


e, (k-1) = u(k-1) - y(k-1) (7.4a) 
e,(k-1) = e, (k-1) -37 e} (k-1) (7. 4b) 
yk) = .02e,(k-1) - .01889e,(k-2) + 2.0y(k-1) 

- y (k-2) (7.4c) 


Manipulating equations (7.4a)-(7.4c) gives 


Wik) = .02u(k-1) -.01889u(k-2) +1.98y(k-1) - .98111y (k-2) 
Ree Cen (eo oosi4e333u> (ko) 
Pe O266ey > (k-1)) = .003148333y~ (k=2) Ges 
2 2 
ey, One 1 ) nae ogaaisu (kaa) (e=2) 


Ramen thy “eed 0944 5a (k= 2) y > (k=2) - 


Equation (7.5) indicates which non-linear terms of 
equation (7.1) should be retained. Therefore, the third 
order non-linear approximation model should contain the 


following terms: 


Meat, U(k-2), y(k-1), y(k-2) linear terms 7G) 
3 3 E 
Mete—1), u (k-2) input cubic terms C7365 } 

3 3 , 
Mme 1l), y (k-2) SutouE cubic terms” (7.6c) 


Ti 





p 2 
@ (k-l)y(k-1), u (k-2y(k-2) third order cross 

72 6d,) 
PEoaduGtmeerins 


Metk-1)u(k-1), y*(k-2)u(k-2) 


The AKI algorithm is primarily modified to include the 
hybrid signals (7.6) in the structure of H(k), the measure- 


ment vector. Since now the measurement vector, H(k) takes 


the form, 
Atk) = (u(k) ... u(k-m),-y(k-l) ... -y(ken) ,u>(k) 
: u> (k-m) ,-y> (K-1) eo ~y? (ken) ,-(u (k-L) y(k-1))... 
2 2 2 
Po Oem cem)) (u(ket)y  (kel)) ... (u(ken)y (k-n))] 


ey) 


the AKI algorithm calculates the coefficients associated with 
the special case of the generalized non-linear ARMA model, 
Seuation (7.1). 

Romeo using a similar approach [Ref. 7] computes a least 
Meares Curve fit for the third order truncation of the 
Taylor series for the sine function over the interval 


(0,7/2). Beginning with the approximation, 


Sine) | = Sees ooneree (7.8) 


Romeo finds the values for a and 8 to be .865 and -.095 
respectively. Performing the same operations as those used 


in obtaining equation (7.5) we arrive at, 
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ic) = .020u(k-1)-.01889qu(k-2) + (2.0 - 020) y(k-1) 
+ (.018890 -1)y(k-2) - .028u7(k-1) + .018898u° (k-2) 
23 2 23 
= .O1889Ry" (k-2) + .068u-(k-l)y(k-1) + .028y7 (k-1) 
~ .056670u~ (k-2) y(k-2) - .068y*(k=1)u(k=1) 


+ .0566708y~ (k-2)u(k-2) (7.9) 


Table 7.1 tabulates the resulting weighting coefficients 
using the three methods just described: (1) analytically 
calculated values using the truncated Taylor series, (2) least 
Squares estimate of the third order sine approximation, as 
computed by Romeo [Ref. 7] and (3) the coefficient estimates 
as computed using the AKI algorithm. 

The tabulation clearly shows that the AKI outperforms 
the first two methods overall. That is, the linear terms 
were identified without question; however, the AKI failed 
to identify the coefficients associated with the Eee 
we k-2), mee) y(k=-2) and Go (OD u(K=-2) terms. The reason 


for the failure is not Known and was not investigated. 


ripe, 





Table 7.1 


VON-Cilbahewe i GhueniG  CORPFICIENTS "CALCULATED 


CEENG THREE METHODS 


(1) (2) 
TAYLOR SERIES LEAST SQUARES EST. 


a eS Silat) = ox + 6x 
fee) ~ X— 3X a = .865, 8 = -.095 
y(k-1) 1.98 1.982700 
2) -.98111 ~.98366015 
u(k) 0.0 0.0 
u(k-1) 0.2 .01730000 
u(k-2) ~.01889 ~ 01633985 
uw? (k) 0.0 0.0 
A k-1) = Ogee 00190 
BE (k-2) .003148333 -.00179455 
fe k-1) 003333  -.00190 
y? (k-2) ~.003148333  .00179455 
u (k-L)y(k-1) —.01 - 005700 
wu (k-2)y(k-2) -.09445 00538365 
yi (k-L)y(k-1) -.01 005700 
y“(k-2)u(k-2)  .09445 ~ 00538365 
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(3) 


AKT USING ACTUAL 
Sin(x) data base 
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003278 
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VIM. FINDINGS AND CONCLUSION 


The performance of the Adaptive Kalman Identifier was 
compared to the LMS Adaptive Filter and the Adaptive Recur- 
Sive LMS Filter using data derived from the models discussed 
in Chapter V. Results for the case where the order of the 
model (that is, m and n of equation (2.1b)) are assumed 
known are presented first followed by the analysis of the 
PLL data. The findings for the overmodeled case are presented 
next. The graphs show typical runs and do not represent 


ensemble averages. 


me ORDER OF THE UNKNOWN SYSTEM IS KNOWN 
ie AKI vs. LMS Adaptive Filter 
Using the MA form of the Adaptive Kalman Identifier, 
its performance can be compared against that of the LMS 
adaptive filter. Synthetic data derived from a plant whose 


meadnster function is, 


1 Z 


Mieco + 2.02 °° + 3 0z (8.1) 


was used. However, a fair comparison necessitated that the 
LMS Adaptive filter be "tuned" by adjustment of the conver- 
Semce factor, K 5: Several values for K in the range [-.600, 
-.200] were used. The objective of tuning the LMS filter was 
to achieve a fast convergence time with little or no steady 
State error while not compromising filter stability. The 


three filter weights were normalized and plotted for five 


jal 





convergence factor values. It can be seen from Figures 8.1- 
8.3 that convergence is essentially reached by step thirty- 
five, (k ~ 35). As the convergence factor is further in- 
creased it can be noted in Figure 8.4 that the filter weights 
become more noisy, and that when K. = =.600, filter stability 
is being compromised (Figure 8.5a). 

Using the same data, the AKI converges to the MA coeffi- 
cients in less than five lterations. Referring to Figure 
8.5b, 1t was also noted that as the measurement noise v(k) 
was increased, the LMS algorithm yielded more noisy estimates 
whereas the AKI tended to compensate for measurement noise. 
This iS intuitively reasonable since the AKI incorporates 
into its algorithm the effects of measurement error due to 


noisy sensors. Comparing Figure 8.6 with Figure 8.7 the 


latter figure clearly shows the faster convergence of the AKI 
to the plant coefficient. 

It was shown in Chapter III that the AKI gains would 
approach a steady state value and indeed they did, as Figure 
8.8 shows. However, the instantaneous Kalman gains display 
a more erratic pattern as shown in Figure 8.9. The informa- 
tion of Figure 8.8 was gleaned from Figure 8.9 by computing 
the average of the individual gains at each iteration using 


the algorithm, 
E{G(k+l)} = E{G(k)} + goyIG(k) -E{G(k) 1] (8.2) 
where, 


E{G(k+1)} is the one step prediction of the 
average value for the gains, 
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Figure 8.7. LMS Adaptive Filter Weights 
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AKI Moving Average Instantaneous Gains 





E{G(k)} is the present average value for the gain 


G (kK) is the present value of the instantaneous 
gains 

+ is the averagin ain 

eel ging g 


This method of obtaining a "running" average is generally 
well known, (see for example Ref. 39). The example presented 
was not the only one used but a good representative of the 
operation of the AKI when the data was derived from a moving 
average process. 

@- AKI Applied to Autoregressive Data 

Before applying the AKI to identifying the coefficients 

of a complex ARMA plant, it was first tested using data 


derived from a plant whose transfer function specified an 





autoregresSive process. The transfer function used was, 
IL 
Goo SS Re a 
26.0.=1 9.0 -2 iO =3 
BeO Tae Oe 1 24) 0" 
a 
0 
ee (33) 
-1 -2 -3 
1 - b,2 ~ Eee - bz 


The AKI had no problem converging on the plant coefficients 
(Figure 8.10) producing the following estimates at the yee 
Meetation (Table 8.1). It essentially converged on the plant 
coefficients in approximately five iterations. The gains 
computed by the AKI for this case also displayed convergence 


to a steady state value. The gain history for this example 


ms depicted in Figure 8.ll. 
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Figure 8.10. AKI Autoregressive Plant Coefficients 





Figure 8.11. AKI Autoregressive Averaged Gains 
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BCTUAL VS AKI ESTIMATES 


ACTUAL 
Ao 1.0 
by =r 08333 
b. =. 3750 
b. -0.04166 


Table 8.1 


POnmCOunmlGlONio OF EQUATION 8.3 


A ESTIMATES 
Or2995 

= 10e5 

=O) 0 


SOS OMS IL 


3. AKI Applied to ARMA Data 


The next logical 


step was to use the AKI to identify 


the coefficients of a general ARMA process. One of Perry's 


models [{Ref. 6] was used 


HOmmenNas DUbDOSe., Specifically, the 


transfer function of the plant was, 


Miz) = 


i 


eos Oe 2 - b. 


a 


a) eee 
fe 1d4c8 2 1 .454Gn 


- - 982° 
a= Jeeaso2 ~ +.407452 — 
een age 
1 2 
(8.4) 
—-.7 - b aoe 
Ss 4 


The ARMA plant was subjected to the same conditions which 


Perry [Ref. 6] describes. 


That iS, unit variance, zero mean, 


white gaussian noise was used as the input. The reader is 


reminded that the output signal processed by the AKI was cor- 


rupted by measurement noise as described in Chapter V. The 


Output data used in Perry's examples, however, reflects noise- 


less measurements. Table 8.2 tabulates the results for the 


coefficient estimates computed by the AKI. Even though the 


results presented in Table 8.2 represent the coefficient 


oJ, 





Mable 3.2 


PenuaAt VS AKI ESTIMATES FOR COEFFICIENTS OF EQUATION 
Seo oe bol frERATION 


ACTUAL AKI ESTIMATES 
Ay i. COCO Osi (0) 
ay 1.40000 1.40100 
a. 0.98000 0279260 
a 0.00000 Oy012 62 
a, 0.00000 = OMG0922 
by 1.14000 1.14000 
bd. alg Bee Sk ASE UN 
b, 0.88490 .88740 
Dy ~-0.40745 = 202 30 


estimates at the ayia iteration, 1t can be seen from Figure 


8.12 that the AKI has essentially converged by the Maes 


lteration. We note also the characteristic convergence of 
the averaged gains to steady state values in Figure 8.13. 


As a means of comparison with Perry's results, the poles and 


h fe 


zeros of the AKI estimates at the ee , 90 and oe itera- 


ti0n are plotted in Figure 8.14a and Figure 8.14b. Perry's 


results Figure 3.7 [Ref. 107,208) for his lattice 


6700. 
modeling of the plant represented by equation (8.4) are 
reproduced for convenience. 

As was noted for the previous cases, the instantaneous 
gains appeared erratic, Figure 8.15, whereas the averaged 


gains converged to some steady state value. 
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Figure 8.15. Instantaneous Gain Values, AKI(4,5) Model 


4. AKI vs Adaptive Recursive LMS Filter 
Using Feintuch's proposed algorithm [Ref. 18] and 
repeating the simulation presented in the rebuttal to Johnson 
and Larimore [Ref. 19], a comparison was made between the 
Operation of the AKI and the Adaptive Recursive LMS filter. 
The Adaptive Recursive LMS filter required "tuning" of the 


convergence factors, Ke and k to what appeared to be the 


b’ 


Optimal performance features of (1) fast convergence and 
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(2) stable estimates. The "tuning" process was a trial and 
error procedure beginning with the convergence constants 


given by Feintuch, 


Bans ig (8.5a) 


n~ 
Il 


k fer. 100 (8.5b) 


b 


The Adaptive Recursive LMS algorithm was then implemented 
to estimate the coefficients of the plant whose transfer 


function was, 


-l 
| 0.05 - 0.40277 ais 
1 lo Bc le Het 62 Se 
a vi 
See (S06) 
lect —ae 2 ne 
Jk 2 


The best response obtained by trial and error resulted in 


mBme Convergence constants, 


i ee ane (ES 


r 21s), 2) Salon ae Ca 


b 


Feintuch reported that after Spaloe iterations the estimates 
had converged on the coefficients of (8.6a) with a .2096 
normalized rms error. Using the convergence constants (8./7) 


convergence was essentially reached by the 4,000°%% iteration. 


Feintuch reports the resulting estimates for the coeffi- 
cients at several intermediate iterations from 8,192 to 65,536, 
however 8,192 was the minimum number of iterations reported. 


og 





It was felt that the operation of the adequately tuned 
Adaptive Recursive LMS filter could be fairly compared with 
the operation of the AKI. 

Table 8.3 tabulates the performance of the two system 
Mmeentification methods and Figures 8.16 and 8.17 graphically 
present the responses. Figure 8.18 not only shows at what 


point the gains of the AKI reach a constant value but also 


Table 8.3 


ACTUAL VS AKI. AND ADAPTIVE RECURSIVE LMS. FILTER 
Pome wh osORsCOBLE LCIENTS OF EQUATION 8.6 


ACTUAL AKI ESTIMATES (k = 150) ADAPTIVE RECURSIVE 
LMS (k = 3,990) 
Ao O05 Cro 2107 0.05045 
ay -0. 40 -0.4007 10/48) 5) SiG) 
by 1314 ees U0 eZ oO 
b, “0.25 -Q0.25090 -0.2442 


provides a measure of confidence that the unknown system has 
been identified. It is evident from Figure 8.19 that both 
methods seem to identify the poles and zeros of the actual 
Plant; however, the poles computed by the AKI are closer. 
The zeros computed by the AKI are approximately .080 units 
farther from the true zero than is the Adaptive Recursive 
LMS filter estimate. All aspects considered, the AKI takes 


considerably fewer iterations to arrive at its estimate. 
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Figure 8.19. Pole/Zero Models Produced by AKI and Adapt. 
_ Rec. LMS Algorithms for a Plant with the 
Characteristic Transfer Function of Eqn. (8.6) 
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o>. AKI Applied to PLL Data (Linear Region) 


To identify the coefficients associated with the 
general autoregressive-moving average representation for the 
phase locked loop, the input noise signal power (mean square 
value) was kept at (10 deg) *. In this manner, the linear 
region of the PLL was invoked. The input/output data was 
analyzed by the AKI resulting in the estimates shown in 


Table 8.4. Figures 8.20 and 8.21 show the responwof the AKI 


Table 8.4 


ACTUAL VS AKI ESTIMATES OF THE ARMA 
REPRESENTATION FOR THE PLL (LINEAR REGION) 


ACTUAL AKI ESTIMATES (k = 82) AKI ESTIMATES (k = 350) 
Ay 0.000000 -0.0001046 —O 000692 
ay 0.020000 0.0198500 0 .0199800 
a, =). 08390 =O olor 20 -0.0188900 
Dy 1.980000 1.980000 1.978000 
b. -0.981110 -0 .981000 =029790 


when applied to the identification of the PLL data. It was 
noted that though the AKI correctly identified the coeffi- 
cients of the linear PLL, the AKI gains were large due to 
the weak signals (input/output data) incorporated in the 
measurement vector, H(k). 
Beeeeakl Applied to PLL Data (Non-linear Region) 

When analyzing any non-linear system the engineer 

must bring to bear all his analysis techniques on the problem. 


The PLL was therefore studied using classical root locus 
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Figure 8.21. Averaged Gain Values, AKI(2,3) PLL Model 
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techniques to preview the possible outcomes when being modeled 
by the AKI. Stability analysis of a more complex PLL system 
using the root locus technique has been previously presented 
in the literature [Ref. 40]. 

The root locus technique can be applied to the PLL 
presented in this thesis by first assuming that the sin(:) 
mock Of Figure 5.7 is a variable gain, AX. This is not a 
restrictive assumption since the PLL during operation generally 
tracks small deviations. The loop gain can be written by 


inspection as, 


-l, -l 
biz) = AC.02)(1 = .944527") 2° (8.8) 
(ple) acne) 
Wino ze, — .9445)2 
‘pei 


From equation (8.8) the root locus of L(z) is drawn as shown 
in Figure 8.22. Even though the root locus technique is 
generally used when the signals in the system are considered 
Mecerministic, it is not surprising that some of the results 
obtained are nevertheless valid. When a moderately strong 
input noise identification signal [E{u*(u) } AS: deg) “1 was 
used, the pole-zero locations of the PLL seem to follow the 
classical root locus analysis. However, when the input noise 
identification signal power is increased beyond (25 deg) ° 
the pole-zero locations do not follow the expected behavior 
predicted using the root locus method as can be seen in 


Begure 8.22. 
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Me@uge §.22. PLL Root Locus Analysis vs Computer Roots 
of the AKI ARMA Model 
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The departure of the pole-zero behavior from what was 
expected was analyzed by closer investigation of the linear 


terms of equation (7.8). 


y(k) = mUIzZouck—i) =~ .018e9eu(k=-2) + (2.0 —- .02a)y(k-1) 
an 
+ (.01889a -1l)y(k-2) (8.9) 


Recalling that a third order Taylor series approximation of 


the sine is the functional expressed by equation (7.7), 


sin(x) 2 ax + 6x? eg) 


one notes that the linear region is described when a = 1 and 
B= 0. It is therefore reasonable to study the variation of 
a with respect to input noise power. 

The average value a of a, was computed by equating 
the estimated AKI coefficients to the coefficients of like 
terms in equation (8.9) for several input noise power levels. 


The relation used was 


~m nm 


2.0-b 1.0 +b a 


; il 2 1 Cakes 


= [a ¢§ ae tt a = 


. ~02 ~01889 noe OTS89! 


meethe 350th iteration. The relationship between a and the 
input noise power level is readily apparent from Figure 8.23. 
This result is plausible Since if one considers the input- 
Output relationship of the sine block of Figure 5.7 one obtains 
Figure 8.24. Superimposing the gausSian probability functions 


Of the different input noise signals, it can be seen that for 
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low power levels of the input signal, identification of the 
linear parameters of the overall system (equation 5.9) can 
be made. As the input power is increased beyond (25 deg) * 
the input-output characteristics of the sin x block are no 
longer approximately linear, showing its effect in Figure 8.24 
memaedeparture from its linear operation, «4 = 1. Therefore, 
by monitoring a one can determine when the overall system is 
entering its non-linear operating regime. 

Since the same functional dependence between each of 
mae AKI estimates and a did not exist for all of the input 
noise powers considered, knowing a did not provide any infor- 


mation of the pole locations but did provide a measure of the 


degree of non-linear operation. 


B. ORDER OF THE UNKNOWN SYSTEM IS NOT KNOWN (OVERMODELING) 
Meenwl VS IMS Adaptive Filter 
Using the data derived by operating a plant with the 
transfer characteristic of equation (8.1) the operation of 
the AKI and the LMS adaptive filter was compared when the 
orders used in the identifier and the LMS filter were greater 
than the known process that generated it. The following 
overmodeling cases were studied: 
(1) MA model is greater than plant MA process 
(2) ARMA mcdel is fitted to MA plant 
When the MA model order 1S greater than that of the 
plant MA process, it was found that both the AKI and the LMS 


filter would compute coefficients close to zero for the higher 
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order coefficients. This effect can be seen in Figures 8.25a 
and 8.25b when a fourth order MA process 1S used to model the 
actual second order process represented by equation (8.1). 
Since these are a "zeros only" plant and model the over- 
modeling essentially causes zeros to appear at the origin of 
the Z-plane of the model transfer function. Table 8.5 sum- 
marizes the results for two overmodeling conditions at the 250th 
iteration using the AKI. When a purely moving average process 
represented by equation (8.1) was modeled as an ARMA process, 
one must direct his attention to the poles and zeros of the 
model transfer function which the AKI computed and compare 
them to the actual plant poles and zeros. It was not readily 
apparent from the resulting coefficients that the plant had 
been identified. The following example will help clarify 
what is happening. 

For the data produced by the second order plant 
(equation 8.1), an autoregressive moving average (ARMA) model 
of orders 2 and 5 respectively was fitted. It can be seen 
from Table 8.6 that no firm conclusions can be drawn about which 
coefficients actually identify the plant. We, the analysts, 
knowing the form of the plant which produced the data, could 
l’ 5! Ax, ay and ac are small 


enough to be ignored. Hence, we can identify the plant cor- 


qualitatively state that b b 


rectly. However, inspection of the poles and zeros of the 


transfer function which the AKI computed, 
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Figure 8.25a. AKI(0,4) Model of a Third Order MA Process 
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Table 8.5 


OVERMODELING OF A MA PROCESS USING THE AKI (k = 250) 


Actual Model 


Coeffs 
1 
2.0 
S130, 
0.0 
0.0 


01618 


zeros: -1.0+1.414)5 


a0) ce ay 


vy 


3rd Order AKI 


5th Order AKI 


MA Model MA Model 
.9995 9991 
2.0000 2.000 
3.0000 3.000 
.893 x10> 0003191 
= 001082 
— 0003198 
-1.000 +1. 4145 -1.001 +1.4143 
2977 1007 .0227 + .04285 


Table 8.6 
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OVERMODELING OF A MA PROCESS USING THE AKI OF ORDERS 


Actual Plant Ccefficients 


OF 


0 


0 


a0 


AR = 2, MA = 6 


.00467 
SLO) S 010 
1.0000 
1.9940 
2.8840 
=Oee2o7 


—O s3i42 


AKI (2,6) Model Ccefficients 


®.0001517 


IES 





Momo + 2.66402 * = 0.22572 > 
ot: = 0.314227" + 0.015172 ° - 
DeeeO0467izem— —- . 10502 
(Sree) 
results in: 
pedes:  sgex8217 
.3264 
momese., — 999 6mte1 . 4145 
e270 
~.3212 
2B 


The observation to be made is that there are two pole-zero 
combinations which are near cancellation. This suggests that 
the AKI algorithm be rerun with the AR and MA orders reduced 
Dy at least two. 

The implication of the analysis of this section is 
that a model can in theory be found for a given set of input/ 
Output data. Further, a parsimonious model can be identified 
by careful observation of pole-zero combinations which are 
near cancellation and of stray zeros near the origin. 

Z AKI Applied to AR and ARMA Data 

Essentially the same characteristic results found in 
Section VIII.B.1 were confirmed when the data produced by 
plants defined by the transfer functions of equations (8.3) 
and (8.4) were analyzed by the overmodeled AKI. That is, 


pole-zero pairs near cancellation and zeros near the origin 
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were produced by the AKI. Figures 8.26 through &§.28 have 
been chosen as representative pole-zero plots of the transfer 
functions computed by the AKI for the AR plant of Section 


Mert. A.2 and for Perry's model (Section VIII.A.3). 


= CONCLUSION 

This work indicates that the Kalman filter algorithm 
heretofore generally used as a state estimator, or in augmented 
form to estimate parameters (in which case the parameters 
are treated as states), can also be formulated in an adaptive 
manner to iteratively estimate the coefficients of an ARMA 
equation explicitly. This approach, termed the Adaptive 
Kalman Identifier (AKI), summarily identifies the unknown 
System whose input/output data is being processed. The LMS 
adaptive algorithm of Widrow, and its modification by Griffiths 
(in which the convergence factor is selected to be inversely 
proportional to the input signal power) are shown to be sub- 
optimal cases of the AKI. An additional insight provided by 
the AKI is that it indicates clearly how measurement noise 
might be taken into account in the LMS adaptive formulation. 

The operation of the AKI was checked by way of simulation 
and compared with two existing identification techniques: 
(1) the LMS Adaptive nonrecursive (MA) algorithm and (2) the 
Adaptive Recursive ARMA LMS algorithm. It was found that not 
Only are the two LMS filtering techniques special suboptimal 
cases of the AKI; but, further, the AKI exhibits superior 


convergence and modeling properties for the cases where (1) the 
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Figure 8.26. AKI(3,1) Overmodeled AKI (4,2) 
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order of the unknown plant is known and (2) the order of the 
plant 1s overmodeled. Additionally, the simulations indicate 
that accuracies similar to those obtained using lattice 
modeling methods, can be achieved using the AKI at a decidedly 
smaller number of iterations. 

By making minor modifications to the measurement vector, 
H(k), that is by using hybrid signals, the AKI was used to 
identify the linear and non-linear ARMA representations of a 
phase locked loop with success. Interestingly, the AKI tech- 
nique appears to enable one to discern when a potential non- 
linear system enters its non-linear mode of operation, by 
closely monitoring the coefficients of the linear portion of 


the generalized non-linear ARMA model. 


PeetoricsS FOR FURTHER CONSIDERATION 

Several areas for further study directly and indirectly 
related to the AKI were uncovered. Foremost, a rigorous 
convergence proof is desirable. Although the connection was 
made between the AKI algorithm and the initial equation from 
which the lattice modeling algorithm is developed, similar 
comparisons (such as Chapters III and IV) could provide more 
insight into the operation of both. The multichannel AKI is 
the logical development of the single channel AKI presented. 
And, lastly, refinement of the AKI software (using the NL = 5 
option) to include ARMA modeling of time series using a Box- 
Jenkins approach [Ref. 26] is feasible. A limited number of 


Simulations using Monterey rain data and series C [Ref. 26] 


es 





ea seem to indicate that an application exists for the AKI 


mr this area. 
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APPENDIX A 


tooo foc ste Wo LNER PROBLEM 


The discrete Weiner problem considers optimally filtering 
a desired signal from unwanted stationary noise. The cri- 
terion of optimality used is minimization of the mean squared 
error between the output and the desired signal. Generally 
one desires that the optimal filter, which is the device 
being sought, be time invariant and that it be able to accept 
a Signal, s(k), and noise, n(k), wnere each are samples from 
Stationary random processes. In other words, we want a 


Gevice which can accept 


Swe SK) Font) (A.1) 


as an input and produce at its output s(k+A) or some linear 
function thereof where 4 1s some known delay. 

Assuming that the desired output, d(k), is the response 
to a specified sampled data linear system whose transfer 
Bunetion, H4(2), 1s given, then the error between the desired 
Siiepuc and the output of the filter, y(k), we seek is, 


ee) = Wd(k)r = y (kK) ; (A.2) 


Figure A.l depicts the discrete Weiner problem formulation. 
The derivation from here follows the one presented by Maybeck 
(Ref. 27] for the continuous case. From the block diagram 


we have, 
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s (k) 





Figure A.1l. The Discrete Weiner Problem 


CO 


ys) a) Even) (a) (A.3) 


,=—-0 


or equivalently, 


co 


aS) ) x(k-i) h, (i) (A.4) 


L,=-—-0 


Substituting equation A.4 into A.2 and squaring we have, 


Cc 


(ARES ) 


e*(k) = a%(k) = 2a(k) J hy (idx (k-4) 
j=-o 
+ £ ) x(i)hg(k-i)]{ ) x(i)h,(k-i)] 
La ae 


IeZz0! 





Taking the expected value of A.5 we obtain, 


Ele*(k)} = Yagi) - 2 2 bh, (i)¥g, (4) 
1=-0 
(A.6) 
+ 2 ne (i) i stil aca le 


where notation Vay (™) denotes the expected value of the 


mmeauct Of u(k) and v(k+m). That is, 


Yay ™) = Ef{u(k)v(k+m)} = Efu(k)v(k-m) } (Caaes 


otherwise known as the autocorrelation of u(k) and v(k). 
Using variational techniques and letting, 


h, (1) = Drage ela (a) (A.8) 


substituting equation A.8 into equation A.6, taking the 
partial derivative with respect to epsilon, ¢€, and setting 
the partial derivative equal to zero we have: 


foo) 


StL b Bong Ct) Vy (Bod) 7 Vig in) | = 0 (A.9) 
The term within the brackets is the discrete form of the 
Weiner-Hopf equation and must be equal to zero if equation A.9 
1s to be valid since Ah(n) > 0 by definition. Therefore, 


oO 


Pheu, (ni) = Ug im) (A.10) 


a —o 
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Equation A.10 is the most often encountered in truncated form 
in linear prediction theory, Pade approximation, adaptive 
filtering and lattice filtering. The truncated version of 
A.10 which is generally used in solving the discrete Weiner 


problem is, 


Sopot 2) Woe as = Ug ln) (A.11) 


| ~1%2 


1=0 


ein matrix form, 


xx (0) Yg (FD) cee Vy lom) T fag, (0) bg (0) 
Vg!) Vy (0) vee Hy (en Pho) Veg (1) 

=| (A.12) 
Mee) (n=)... (0) Dope (M) Yq (R) 
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